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Abstract. In this paper we consider the Long- Wavelength Approximation Scheme 
(LWAS) in the framework of the covariant fluid approach to general relativistic 
dynamics, specializing to the particular case of irrotational dust matter. We discuss 
the dynamics of these models during the approach to any spacelike singularity where 
a BKL-type evolution is expected, studying the validity of this approximation scheme 
and the role of the magnetic part of the Weyl tensor, Hab- Our analytic results 
, confirm a previous numerical analysis: it is Hah that destroys the pure Kasner-like 

ly-^ ' approach to the singularity and eventually produces the bounce to another Kasner 

. phase. Expanding regions evolve as separate universes where inhomogeneities and 

\ anisotropics die away. 
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bj;). 1. Introduction 

■ Approximations methods are and have been very important in general relativity and its 
' applications to cosmology and relativistic astrophysics. In the case of cosmology linear 
perturbation theory ^12] is at the heart of the present developments, in particular, in the 
analysis of the observed anisotropics of the Cosmic Microwave Background (see, e.g., jHj). 
Another approximation technique widely used in cosmology is the so-called Long Wave- 
Length Approximation Scheme (LWAS). Broadly speaking, it consists in neglecting 
spatial gradients of the variables describing the cosmological models, assuming they are 
small in comparison with their time derivatives. Since the time-scale of variation in 
cosmology is given by the local Hubble expansion rate, the approximation consists in 
neglecting inhomogeneities varying over a scale smaller than the local Hubble horizon. 
One can then try to get information about scales smaller that the Hubble radius by 
computing their effects through a series expansion in spatial gradients, which constitutes 
the LWAS (sometimes also called the gradient expansion method or long wavelength 
iteration). This scheme is then suitable to study large-scale structure formation and 
issues related with it. One important advantage of using the LWAS is the fact that 
neglecting spatial gradients we are still taking into account the non-linear character of 
the dynamical equations. Actually, the dynamical evolution that comes out will be fully 
non-linear and essentially the same as the one for homogeneous cosmological models. 
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In the literature there have been several approaches to the LWAS: using this 
approximation Belinski, Khalatnikov, and Lifshitz (BKL) studied jl] (see also the 
general behaviour of cosmological models near the initial singularity. In [HI Ej it was 
used to study several important issues in structure formation, from the gravitational 
instability mechanism to the evolution of the primordial inhomogeneities. Higher-order 
expansions in the LWAS and their applications have been developed in [H] by using a 
Hamiltonian approach. A closely related iterative method has been worked out in PlllUj. 
Another approximation scheme close to the LWAS was presented in [TT] . 

In this paper we study a different formulation of the LWAS based on the covariant 
fluid dynamics approach to cosmology (see, e.g., [121 Cni EI)- In this sense our work is 
complementary to the perturbative gauge invariant formalism in ^Sj (see also fl] and 
references therein). We focus on the study of irrotational dust models (IDM) with a 
cosmological constant A > 0, and more specifically, on the dynamical aspects near the 
singularities. One of the main conclusions of this study is that the LWAS breaks down 
due to a gravitomagnetic effect. This points at what it should be expected in the light 
of the standard BKL picture: that in the generic case the approach to a singularity 
has a first phase in which the solution goes into a Kasner phase, after which there 
is a "jump" into a transient period in which the solution moves away from a Kasner 
solution, and then goes into a new Kasner phase. With our language, and within the 
LWAS, this jumping is caused by a gravitomagnetic effect that is necessarily associated 
with non degenerate cases (here degenerate will mean two equal shear eigenvalues). 
From the point of view of the evolution equations, the LWAS reduces the quasi-linear 
Partial Differential Equations (PDEs) system of Einstein's equations, or equivalently the 
covariant fluid equations, to Ordinary Differential Equations (ODEs). The evolution in 
this phase proceeds in a "silent" fashion E] approaching a Kasner singularity, till 
the effects of the gravito-magnetic Weyl tensor Hat are no longer negligible and a bounce 
to a new Kasner phase is produced. 

The plan of this paper is the following. In section |21 we review the fluid dynamical 
approach to IDMs and choose the most suitable form of the equation for the application 
of the LWAS. In section El we show how to implement the LWAS within the chosen 
approach and we present the flrst-order solution. In section HI we discuss the validity of 
the approximation from the point of view of the covariant fluid dynamics approach. We 
end with some remarks and conclusions in section El 

Throughout this paper we will use units in which SttG = c = 1. Tensorial 
components will be expressed both in coordinate charts and in vector basis. Then, 
our conventions for indices are: spacetime coordinate indices are denoted by the lower- 
case Latin letters a, ...,/ = 0, 1, 2, 3, and spacetime indices associated with an arbitrary 
basis {eo, ei, 62, 63} by the rest of lower-case Latin letters m, z = 0,1,2,3. Indices 
with respect to an orthonormal triad of spacelike vectors {61,62,63} are denoted by 
lower-case Greek letters a,..., A = 1,2, 3 and spatial coordinate indices by the rest of 
lower-case Greek letters 11, . . . ,uj = 1,2,3. 
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2. Irrotational Dust Models 

IDMs are important both for the description of the late universe as well as for the 
study of the gravitational instability mechanism. Their energy-momentum distribution 
is completely described by the fluid velocity u, a unit timelike vector field = — 1), 

and its associated energy density p. The energy-momentum tensor is given by 

Tab = pUaUb , (1) 

Since there is no pressure, the energy-momentum conservation equations tell us that 
the fluid worldlines are geodesies, i.e., u'^VbU"' = 0. Assuming the fluid flow to be 
irrotational {u[a'VbUc] = 0) implies that the fluid velocity is the normal to a foliation of 
spacelike hyp ersurf aces. Then, there exists a time function r(x") such that u is given 
by 

U = , U = —dT . 

or 

That is, r is at the same time the proper time of the fluid elements and a label of 
the hypersurfaces orthogonal to u, which we will denote by S(ti) = {r(x") = ri : 
constant} . A system of coordinates for these hypersurfaces can be formed from any three 
independent first integrals of u, {y'^{x°')} with u^daV'^ = 0. Then, {T,y^} constitutes 
a set of comoving synchronous geodesic normal coordinates, in which the line element 
takes the following form 

ds^ = ~dr' + h,,{r,y'')dy^dy'', (2) 

where h^^i, are the nonzero components of the induced positive-definite metric (first 
fundamental form) on the hypersurfaces S(r). This coordinate system is fixed up to the 
transformations: r t' = t + constant, and — > y'^' = f^' {y"), being f^' arbitrary 
functions. 

A very convenient description of these models is provided by the covariant fluid 
approach introduced by Ehlers [12] (see [ini Ej for details). For the purposes of this 
paper it will be useful to write the equations using an orthonormal tetrad adapted to 
the fluid velocity: {eo = 11,61,62,63} with • e„, = rjmn = diag(— 1, 1, 1, 1) . For 
convenience we fix partially the freedom in the choice of the triad {e^,} by requiring it 
to propagate parallely to the fluid flow 

e„« = u^VbCa" = 0. 

This choice also makes the local angular velocity to vanish 

= ^£°''^^ei3 • 65 = [ea/ss = Val3So) ! 

where rjabcd is the spacetime volume element. In this framework, the variables that we 
need to describe IDMs are the following, (i) Metric variables. The components of the 
triad vectors in adapted coordinates {r, ?/^}, Ca^. (ii) Connection variables. The spatial 
commutators, 7"/3a, defined by the commutation relations between the triad vectors: 

[e/3, 6a] = 7"/3Aea , 7"[/3A] = 7"/3A • (3) 
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We can instead use, the following variables introduced by Schiicking, Kundt and Behr 
(see, e.g., [H]): 

which contain the same information as 7"/3a- (hi) Kinematical variables. The expansion 
(= VaW") and the shear tensor of the fluid worldlines. The shear is a symmetric and 
trace-free spatial^ tensor 

(Tab = h(^a''hb)'^VdUc - (0/3)/la6 = VaMfo - (9/3)/la6 • 

(iv) Matter variables. The only non-zero component of the energy-momentum tensor is 
the energy density p = Too. The Ricci tensor is, through Einstein's equations, given by 

Rab = P UaUb + |(P + 4A)fi'afe , (5) 

where A denotes the cosmological constant, (v) Weyl tensor variables. The Weyl tensor 
Cabcd describes the spacetime curvature not determined locally by matter fields. Its 
ten independent components can be divided into two spatial, symmetric and trace-free 
tensors 

Eab = CacbdU'^u'^ , Hab = *CacbdU''u'^ i*Cabcd = ^Vab""-^ Ccdef) , 

which are called the gravito-electric and gravito-magnetic fields respectively. Whereas 
the gravito-electric field produces tidal forces, having a Newtonian analogue (the trace- 
free part of the Hessian of the Newtonian potential), the gravito-magnetic field has no 
Newtonian analogue. 

The equations governing the behaviour of these quantities come from the Ricci 
identities applied to u, the second Bianchi identities, Einstein's equations (jSj) and 
the Gauss-Codazzi equations (see [121 Ej for details). To obtain them we have to 
split the covariant derivative Va into a time derivative along the fluid worldlines, 
A"'"'f)... = u'^VcA"'"'b..., and the induced covariant derivative on the hypersurfaces S(r), 
Df-A^-'b-.- = h^-e' ' ' h-^bhc'^Vd.A'^'" f..., where A°'"'f,:. is any arbitrary tensor field. To 
simplify the equation that will appear in this paper it is convenient to introduce the 
spatial divergence and curl of an arbitrary 2- index symmetric tensor§ Aab (see, e.g., |19j ) 

div {A)a = D^Aab , curl Aab = ScdiaD'Ab)'^ , 

where Sabc {^'^^^^def = 3\h'^°'dh''eh'^^ f) corresponds to the volume element of the 
hypersurfaces S(t). The projection of these definitions onto a triad {e^} gives|| 

cuTl{A)ap = e^^^aidixi - aix\)Apys + sA^p - ^n^jAp^s , 

div(A), = {ds - 'ias)A\ - eJ^np^Axs , 
X Spatial means orthogonal to the fluid velocity. 

§ These definitions are analogous to those for vector fields: div(A) = DaA'^ and curMa = SabcD^A'^ ■ 
II Angled brackets on indices denote the spatially projected, symmetric and tracefree part: = 
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where da = Ca^dyt^ . Finally, for two arbitrary spatial symmetric tensors, Aab and Bab, 
we define the commutator as 

[A, B]ab = 2Ala'Bb]c , [A, B]a = ^SabclA, 5]"^ = EabcA^B^" . 

The resulting set of covariant equations is usually divided into two groups: 
Evolutions equations, giving the rate of change of our quantities along the fluid world- 
lines, and constraint equations, relations containing spatial derivatives only. The explicit 
form of these equations has been given in many places (see, e.g., ^21 ^1 E])- In the 
case of IDMs they can be found in [201 E]- In order to set up a suitable framework for 
the application of the LWAS, here we will look at these equations from a different point 
of view than the usual one. We consider as main variables the set {e^'^, Q, CTa/^}, which 
will be the dynamical variables. Their evolution equations can be written as follows: 

e„^ = -(ie5/ + a,^)e/, (6) 
e = -i02-|a°V/. + |A-fi?, (7) 

&al3 = — ©CTq/S — ^Sal3 ■ (8) 

In these equations, and ^Saf3 denote the scalar curvature and the trace-free part of 
the Ricci tensor of the hypersurfaces S(r) respectively. They must be understood as 
given in terms of 7^^ and their derivatives, through the expression 

where F^^ = ■ (Vg^e^) are the Ricci rotation coefficients associated with the triad 
{cq}, related to -f'j^g by ^a^r^;^ = SaeYxp + ^l^eYaX + ^Xelap ■ The remaining variables, 
{7a/3; ^ai3, H^p], which will be called auxiliary variables, can be found in terms of the 
main ones, either from their definitions or from the constraints. The expressions are: 



lUd,y^=[da,dp\y\ (9) 

P=W- -K + fR, (10) 

Eal3 = — Ci <(J<Jl3>5 + ^Sal3 ■ (11) 

Hap = CUTl{a)af3 ■ (12) 

These quantities are subject to the following constraint: 

div(or)„ - fa^e = . (13) 



The remaining equations involved in this formalism can be derived from those here 
presented. It should be noticed that Einstein's equations enter the scheme above only 
through equation (fTUj) . through the particular way p is related to the geometrical and 
kinematical variables. 
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3. Formulation of the Long Wave-Length Approximation Scheme: the 
first-order solution 

As we have already mentioned the LWAS consists, roughly speaking, in neglecting spatial 
gradients with respect to time derivatives. One can find in the literature several ways of 
implementing the LWAS [HI El E] • In the case of IDMs, the spatial gradients correspond 
to derivatives tangent to an intrinsically characterized foliation, namely, the one formed 
by hypersurfaces orthogonal to the fluid velocity. In this sense, this approximation 
scheme is physically well- motivated for the case of IDMs. Based on this fact, here we 
will adopt a new point of view based on the covariant fluid approach introduced in the 
previous section, which will allow us to clarify several geometrical and physical aspects 
of the LWAS. Although the formalism we develop can be applied in different contexts, 
here we focus in studying the approach to spacelike singularities, complementing in this 
way the work done by Langlois and Deruelle [TU] . 

Let us consider the dynamical equations governing the behaviour of IDMs shown 
in the previous section. In order to apply the LWAS ideas it is very important to realise 
that in the evolution equations for the main quantities, Eqs. (jE))-®, only second-order 
spatial gradients appear, and they are encoded in the spatial curvature tensor, or more 
speciflcally, in and ■ Therefore, the zero- and flrst-order solutions can be found 
in one iteration (we will separate them later) , which means that we can look directly for 
the flrst order solution just by neglecting in these equations terms of order greater than 
one in the spatial gradients. In our case, this means to neglect the spatial curvature: 

^R^O, ^Sab^O, (14) 

which implies, through Eq. ()12|) . that 

DaHbc ~ =^ dw{H)a ^ and cmlHab ~ . (15) 

The flrst important consequence of (fT^ is that from the evolution equations (jH))-® we 
obtain a closed system of equations for B and cTq,^ . 

e = -i02_|c,"V^ + |A, (16) 

O'a^ = -6cr„/3 . (17) 

Second, the expressions for p and Ea/3 become purely local 

p = |02-|a"^cr„;3-A, (18) 

EaP = l^CTaP - (y<J'<yp>X ■ (19) 

In order to integrate equations (fT^ . lfTTj) we will further specialize the triad {ea} by 
taking into account the special structure of the flrst-order equations and the freedom: 
^a' = -^a'^^a with ^q,'" = . We use this freedom to have an eigenbasis of the 
shear. That is, 

cr„^ = for a ^ /3 , (20) 
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This is possible because at the first order we have '^Sa/s = . In [221 (see also it 
was shown that a sufficient condition for this is Hat = 0, and in [25 this result was 
extended to the case in which Hat is transverse, i.e. div{H)a = 0. We have already seen 
that this last condition holds at ffist order [see Eq. (fT3j)]. As a consequence of and 
f|TTHl . Eai3 is also diagonal. Hence [a, E]ab = . 

We can now solve the equations for the ffist-order solution of the LWAS. The 
procedure that we will follow is very close to the one used to find the complete class of 
IDMs with flat spatial geometry [20] • We begin by solving equations (fTH|l . (fTTj) for the 
expansion B and the only two independent components of the shear, whose information 
can be encoded in the following two quantities (see, e.g., and references therein) 

a_ = ^{a22-as3). (21) 

The equations for (0, cr_,_, cr_) then are 



1^11 



e = 



-^(62 + 4 

-ecr±. 



a 



3A) 



and the solution can then be expressed as 



0"_ 



(22) 
(23) 

(24) 



where V = V{T,y^) denotes the proper volume associated with the fluid flow of u and 
E-t are scalars independent of r. From ()22l23p the equation for V is 



2VV -V^- 3AV^ + 3S2 = , 



(25) 



and is a ffist integral of the system 



In order to write the solution of this 
equation one has to consider separately the following cases: (i) A > 0. Here, we have 
to consider also different cases according to a criterion that involves A, S, and an 
integration function that we will call M. Then, when — 4AS2 > 0, the solution of 
()25|) can be written as 



V{t,y^ 



A 



VM2 - 4AE2 sinh 



(r-T) 



sinh 



r - r + ST) 



where T is an arbitrary function of the comoving coordinates y'^. When M'^ — 4AS2 
we have 



V{t,y^ 



M , 



2A(r-r) 



and finally, when M'^ 



V{t,y^ 



-^AAT?-M^ sinh 
A 



1 

4AS2 < 0, the solution is: 



3A 



cosh 



3A 



r - r + ST) 



Where the function (5T, in cases ()26I28|) . is given by 



5T 



In 



M + 2v^S 



M-2VAS 



(26) 
= 

(27) 
(28) 

(29) 
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M 



(ii) A < 0. In this case the solution for the comoving volume is 
V{t, y^") = -^VM2 -4AS2 sin — ^ — (r - 7 ) cos 
where 

01= , arcsm 



T-r + 6T) 



V=3A -4AS2 
(iii) A = 0. The solution for this case is: 

V{t,y^) = lM{r-r){r-r+ST) 

where 

6T 



(30) 



(31) 



4S 



(32) 

M is a time-independent 

(33) 



3H. 



In all these expressions for the comoving volume V{t, y^) 
scalar related to the energy density through equation (fTUj) : 
M 

and a first integral of the energy-momentum conservation equation p + Qp = 0. 
Taking into account that the proper volume (jM]) is defined up to a multiplicative time- 
independent scalar [see Eq. ^^], we naturally choose M to be a positive constant, and 
in fact we could even choose its value. On the other hand, the proper volume is related 
to the local scale factor a(r, y^^) by the expression 

V _ a 
V~a 

where H is the local Hubble parameter. 

Since we will be mainly interested in the behaviour near the singularities we can 
neglect the effect of the cosmological constant and consider only the case (iii) given by 
equations ()31I32|) . This well-known fact can be checked from the expressions we obtain 
for V{t,y^), from where one can see that in the limit r — >■ T all the cases behave like 
the case (iii), i.e., ~ r — T. In any case, the procedure we will follow here to solve the 
different equations can be applied from the beginning to the end to all the cases with a 
non-zero cosmological constant. 

Then, as we can see from fl31l32j) . and assuming without loss of generality S > 
(^ 5T > 0), for r > T we have an expanding fluid (note that this depends on which 
fluid element we are looking at, i.e. on y^), whereas for r < T — 5T we have a collapsing 
one. The region r G (T — (5T, T) is excluded by the requirement of having a positive 
energy density. Then, for an expanding IDM, T{y^) is the Bang time corresponding to 
a fluid element with comoving coordinates y'^. In the same way, T — 6T is the Crunch 
time for a contracting IDM. 

The next step is to solve the equations for the triad vectors {gq,} [Eq. ©], which 
become 



)el^ 



Covariant Fluid Dynamics: a Long Wave-Length Approximation 



9 



e2'' = -|(e + a+ + y3(T_)e2'^, 

e3'^ = -i(e + (T+-y3(T_)e3^. 

The result, in the case A = , is 

= 6„/^(r - r)~P-{T -T + 5T)-^- , (34) 
where the quantities ha^ in fl34p are the time-independent components of a triad defining 



a Riemannian 3-dimensional geometry whose metric tensor is given by q^^p = daph'^pb^v, 
being 6"^ the inverse matrix of ha'^. Moreover, = 2/3 — and 

v^S-2S+ S+ + V3(S + S_) E+ + V3(S-S_) 

3V3S ' = ' = ' ^''^ 

are defined in the analogous way as they are defined for Bianchi I IDMs. Indeed, we 
can check that pa and qa satisfy the Kasner relations 

3 3 3 3 

I] Pa = I] Pa = 1 H ?a = II ?a = 1 , 

SO we will call them local Kasner coefficients. It is convenient, for practical purposes, 
to consider the fact that one can assume, without loss of generality, that in an open 
domain of a given fluid element located at y^, the Kasner coefficients can be chosen to 
lie in the following intervals: 

pig(-|,o), p2e(o,|), p3e(|,i). (36) 

Then, looking at (j34p it is clear that the first-order solution has an initial singularity 
(in general not simultaneous) of the Kasner type, also known as a velocity-dominated 
singularity j2Sl- 

Summarizing, expressions (j23), ()31I32|) . and solve the first-order evolution. 
From them we can find the time dependence of all the variables. With regard to 
the spatial dependence, which is encoded inside the functions T, {6a}, and S-t, the 
only restriction comes from the constraint (fTSj) (the momentum constraint) which is 
equivalent to the following three equations 

9i((T+ + 0) - 3aicr+ - ^3^230-- = , (37) 

d2{(T+ + V3(T_ - 26) - 3a2(a+ + V3(T_) + V?,nvi{V^(T+ - (t_) = , (38) 

a3(a+ - VSa- - 26) - 3a3(a+ - VSa.) - Vlni2{V?,a+ + a_) = , (39) 

To see which kind of restrictions these equations impose we need to compute the 
quantities and Uap, or equivalently the commutators 7"/3a- Using the expressions (jHI) 
we get 



1/3 



'?/3 



(40) 
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where 7^^^ are the (time-independent) commutator functions associated with the triad 
{6q,}, underhned indices do not follow the usual index summation convention, and we 
have used the following definitions: 

d'a = ba^dy^., P = T-T, Q = T-T+5T. 

Then, from ()4p4()|l we have 



<+2 



n- 



P 



d' P 

' P 



Q 



if a = P . 



n 



af3 



{n'^-le^d'^[ip^-pp)ln^]]p-''^Q-'^ (a ^ 5 ^ P) ifa^/3. 

Introducing all these expressions into the constraints ()H7IIH9|1 . we can see that they are 
equivalent to the following three time-independent equations 

d\ [S+ + f M(T - i5T)] = 3a;E+ + VSn'^^E. , (41) 
d'2 [S+ + V3S_ - 3M(r - |5T)] = 3a2(S+ + V3S_) - V3n[.;^{V3J:+ - S_) , (42) 

d'3 [S+ - - 3M{T - l6T)] = 3a^{i:+ - V3S_) + V3n[^{V3^+ + S_) . (43) 

This is not an obvious result since equations (p?fj) - (j!?n|l are linear combinations of 
functions depending on {r, y^} with coefficients depending only on {y^} and, in 
principle, one would expect that from each of them we would get several equations. 
What we have seen is that we only get one for each of them. Equations (pT|) - (jiH|) are 
then constraints on the quantities T, {6^} and T,± . 

At this point, the construction of the first-order solution is essentially finished. As 
we have seen, the time dependence is explicitly known whereas the spatial one is encoded 
in the quantities T, {6a} and T,± . Now, we can split the zero and first order terms by 
taking into account that the zero order is characterized by the absence of any spatial 
gradient. Then, the splitting is given by the relations 

ba^ = 6^^ + ifea'^ , (44) 

S± = oS± + iS±, T = „T + iT, (45) 

where o^± and qT are constants. From we deduce that the commutator functions 
7^% or equivalently, a'^ and n'"'^ do not have a zero-order part. Moreover, it also gives 
a splitting of the spatial derivative d'a 

d'a = dy, + ib^f^dy, . (46) 

Then, it also follows that equations (j4H) - (j43j) are automatically satisfied at zero-order. 
At first order they are 

di + S] = 3a'ioS+ + V3n23oS_ , 

82 + V3iS_ - 2E\ = 3a'2(oS+ + Vs^S^) - V3n[^{V3oJ:+ - oS_) , 
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iS+ — v^iS_ — 2S 



where 



2 



1 



It is worth noting that our zero order solution is exactly Bianchi I, and so is the first 
order part from the point of view of the evolution equations. 

To complete the construction of the first-order solution, we give the expression for 
the remaining quantities, that is, the gravito-electric and -magnetic tensor fields Eap 
and Haji . As we said before, equation (fTn|) tells us that Eap is diagonal since aap is 
diagonal. Then, we only need to compute the quantities and E_ defined as in (j^IJ. 
The result is 



E, 



a 



E. = \{Q- 2a+)a_ = ^ {|M(r - T + \5T) - 2E+} S_ . 

The gravito-magnetic field has five independent components, which can be given 
in terms of the quantities H± [defined as in Eq. (EH)], and = V3H23 , H2 = V3Hi3 , 
and = \^Hi2 . Their form is given by 



frill 
2S 



0-+ - ^(^22 - ri33)0"- 



-n;iP2(Pi-i)g2fe-i) _ 2E_ r ,^^2fe-i)g2fe-i) _ ^/^p2(p3-i)g2fe-i)| ^ 



H 



^{n22 - ^33)0-+ + - 2^22 - 2n33)cr_ 



2S 



2E_ 



p2(p2-l)g2(q2-l) _ ^/^_^p2(p3-l)g2(g3-l)| 



Hi = - V3n23cr+ + {di - ai)cr_ 
= 3|^{9i(iS_)-„S_a;-v^„E_ 



20^- 



'23 o^- y I 2 

(aipi)ln^ + (pi-l)^ + (gi-l)^ 



,S_^ + ^„S+5i[(p2-P3)lnf' 



^^2 = 1^13 (cr+ + V30"-) + 1(^2 - a2)(V30"+ - Cr_) 

= siz {^2(v^iS+ - iS_) - (v^oS+ - oS-)a'2 + 3(oS+ + y3oS-)ri'i 



(v^oS+ - oS-)^ + + v^oS_)92 f(pi - P3) In ^ 



Q 



13 



[(^2^2) In f + (p2 - 1)¥ + (g2 - 1)^] } p-P^-iQ-^^-1 

^^3 = |'^i2(o-+ - V^(y-) - \{dz - a3)(V3cr+ + 0--) 

{-93(V3iS+ + iS_) + (v^oS+ + oS_)a'3 + 3(oS+ - VSoS-)^; 



2 

3M 



(v^„S+ + „S_)^ - f - v^oS^)93 f(pi - P2) In 



Q 



(93j93)lnf + (P3 - 1)^ + - l)^]}p-^'^-lg-'?^-l 



''12 
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Table 1. 


Asymptotic behaviour of the first-order solution. 


Quantity 


Behaviour as t ^ T 


Behaviour as r — > oo 
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In order to understand the dynamical behaviour of IDMs using the first-order 
approximation of the LWAS, we have computed the asymptotic behaviour in time, both 
near the Big-Bang singularity (r —>■ T) and for large future proper time (r — > oo) of the 
main physical quantities that are used in the covariant fluid approach. The results are 
summarized in Table ^ These results represent the generic behaviour of IDMs, there 
are particular cases, characterized by the vanishing of some of the connection quantities 
naf3 or some the components of the gradient of T, in which the behaviour indicated in 
Tabled has to be revised. For the sake of brevity we have not included these particular 
cases here. Moreover, to obtain this asymptotic behaviour we have used the fact that in 
the neighbourhood of a fluid element one can always assume that the Kasner coefficients 
are distributed as in 

4. On the validity of the Long Wavelength approximation scheme 

The LWAS is based on the general idea that spatial gradients can be neglected with 
respect to time derivatives. This assumption puts restrictions to the validity of the 
expressions we have given above. This issue was already studied by Deruelle and Langlois 
in Their criterion to consider the approximation as valid is to compare the third- 
order terms arising from the first-order solution and to check whether or not they are 
smaller than the first-order ones. To that end, they computed the neglected terms from 
the first-order solution, essentially the spatial curvature terms and ^S^^. Then, they 
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realized that in general the effect of the local anisotropy makes some of spatial curvature 
components to blow up. Then, they argue that when this happens one shall recover the 
oscillatory BKL behaviour. 

Our calculations agree with their results. First, for late times (r — oo), the system 
evolves towards an Einstein-de Sitter model as expected (see Table P), whereas for A > 
it would evolve toward the de Sitter model. On the contrary, we have found that indeed 
some components of the spatial curvature blow up in general as one approaches the 
initial singularity (r — T — ^ 0). Nevertheless, here we want to take a complementary 
point of view to deal with this issue. In particular, we only need to analyze the quantities 
used in the covariant fluid approach constructed the first-order solution. In our scheme 
there is a set of main variables, namely {cq'^, 9, cTq/j}, and a set of auxiliary variables, 
{"f^l^, p, Eai3, Ha/s}. The idea being that the secondary variables can be computed from 
the main ones. Of particular interest are the gravito-electric and -magnetic fields. Their 
first-order form has been given in Table ^ As we can see, their behaviour near the 
initial singularity is completely different. While the normalized^ gravito-electric field, 
a dimensionless quantity, is finite there (but depending on the fluid element 
we consider), the normalized gravito-magnetic field, Q~^Hai3, has components that blow 
up. These divergences constitute an indication of the breakdown of the approximation 
scheme as we approach the initial singularity. Here again it is the local anisotropy that 
produces the blowing up of the gravito-magnetic field: this first-order quantity becomes 
dominant over the zero-order ones and the scheme breaks down. This can be seen in 
Tabled where the dependence of the components of Q'^Hap on the Kasner coefficients 
Pa is explicit. In contrast, the gravito-electric field does not depend on them and its 
normalization is regular. Moreover, from Tabled we also see the fact that the effect of 
matter, described by the normalized energy density G^^p (essentially, the cosmological 
density parameter fi), can be neglected near the initial singularity, a well-known result. 

Another interesting difference between the behaviour of the gravito-electric and - 
magnetic tensors is their local/non-local character. Whereas is completely local at 
first-order, that is, for a fixed fluid element ?/" it only depends on the value of the main 
variables at that point [see Eq. (HH)], H^p does depend on spatial gradients of the main 
quantities. More specifically, on the shear-rate spatial gradients [see Equation (|T2|l ]. In 
this sense, and in the spirit of the LWAS, Hab is a purely first-order quantity and the 
spatial Ricci tensor is a purely second-order quantity. The blowing up of H^p (at first 
order) and "^Rap (as computed from the first-order solution) are not independent facts. 
To understand this we have to look at the evolution of the spatial Ricci tensor: 

^R= -lQ'^R-2a^P'^Sap, (47) 

""Sap = |e -'Sap + a^J -'Sp^s - IcTap 'R + cnT\iH)ap . (48) 

In IDMs the gravito-magnetic field has its origin in the anisotropy, Hab = cnTl{a)ab, and 
it is the only source for the Ricci tensor as we can see from equations ()47p48|l . First, 

% Note that we are normalizing using the expansion, which is related to the Hubble length, the natural 
scale of the system, hy Q — 3H . 
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it can generate ^Sab through equation (|iH|) . and ^Sab is a source for the scalar spatial 
curvature in equation (jlTj). 

On the other hand, the dynamical picture we have when approaching the spacelike 
singularity is essentially the BKL picture, in which the system has a Kasner-like 
evolution with some determined Kasner coefficients, say p'^, and at some point the 
system jumps to another Kasner phase characterized by other Kasner coefficients, say 
p'^. Within our approach, we can define the following time-dependent Kasner coefficients 

V3(T - 2cr+ (T+ + V3((T + cr_) ^ _ a+ + y/3{a - a^) 

— ■ l49j 



Pi 



For the first-order solution, these become the time-independent ones defined in fj35|) . 
because all the components of the shear have the same time dependence. To see how 
they change at higher order within the LWAS, it is enough to consider the evolution 
equation for the shear (jHl) with the trace-free part of the three-dimensional Ricci tensor, 
(which is a second-order quantity) computed from the first-order solution. Now, 
in order to write down the evolution equations for the shear it is important to consider 
some important facts. First, that we are considering a basis in which the shear tensor 
diagonalizes, so it can be described only by the scalars ((T+,ct_). Second, in the 
computation of ^Sajj one can see that the diagonal terms dominate over the non-diagonal 
ones as one approaches the initial singularity at r — T. Therefore, for our purposes we 
can consider that '^Sa/s is diagonal and described by the two quantities 



P3 



3c — 3 3c 3c —V3(3q 3c 

0+ = —2 Oil , "J- = —{ 022 — O; 



33) 



This also means that within this approximation the shear and the trace-free Ricci 
tensor commute, which in turn implies that the gravitomagnetic field is transverse: 
dw{H)a = 0. Taking into account all these considerations, the evolution equations 
for the Kasner coefficients ()49|1 can be expressed in a compact way by introducing the 

And the evolution of p is given by 

(50) 



following vector: p = {pi,P2,P3 
p = ^lAp, 



where 



n = 



0"4 



- a_ 3 c 



A 



I 




1 

-1 



1 \ 




1 



-1 







where ^5'-|-+ are computed from the first-order solution, and a± from the zero-order one. 
Equation ()50|) can be solved for each fiuid element and the solution can expressed in the 
following form 



where 



1 



cos 



sm 



v^/(r,y^)l ^Po(y^),(51) 



I{r,y^)^ r n{r' ,y^)dr' , A' 

Jrn 
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From the first order solution we find that in the generic case, i.e. without extra 
assumptions, the quantity /(r, y^) behaves, near the initial singularity, as 

Then, equation ()51|) shows the change that the Kasner coefficients will suffer when the 
spatial curvature is taking into account. However, this equation will not be valid once 
the Kasner coefficient have changed their value significally, since we are approximating 
the spatial curvature from the first-order solution, where the Kasner coefficients are 
static. This means that the behaviour of the spatial curvature along the different 
spatial directions does not change. In order to take into account this fact one should 
compute the spatial curvature from equations ()47p48j) . which describe the evolution of 
and ^Sa/s, mainly driven by a second-order quantity in the LWAS: cm\{H)ab- This in 
agreement with the observation in numerical simulations of the homogeneous Bianchi IX 
models, which exhibit a BKL behaviour, that the gravito-magnetic tensor Hat generates 
the transition between different Kasner epochs (see PT]). 

Finally, we want to remark that the discussion presented here is based on the generic 
solution of the LWAS. One can restrict the free functions and parameters that appear 
in the first-order solution in order to make it well behaved near the initial singularity. 
More specifically, we can impose the following restrictions on the first-order solution: 

T = % : constant , and n'^]^ = . (52) 

The first condition means that the initial singularity must be simultaneous. Actually, 
without loss of generality we can choose % = 0. The second condition in is a 
restriction on the triad {6q,}. Then, under conditions ()52j) one can check that 

^-^0, ^ 0. (53) 



02 ' 02 

With respect to the results shown in TableQthe only thing that changes is the behaviour 
of Hab near the singularity: the precise behaviour can be found in Table El Although 
for the component Q'^H^ the expression in Table El suggests a logarithmic singularity 
for the case ps = 1 , pi = P2 = 0, this is not the case. For this very particular case, 
the coefficient T-l^{y'^) vanishes identically for that situation and hence the result 
is general for the particular case ()52|1 . 



5. Conclusions and Discussion 



In this work we have considered for the first time the application of the LWAS in the 
context of the covariant fluid approach to cosmology [131111. This produces a system of 
ODEs for our main variables {ca^, 9, Ca/s}, giving a "silent" evolution [1^1117] that, for a 
collapsing fluid, proceeds towards a Kasner singularity. However, our analytical results 
show that at a certain point the gravito-magnetic field Hab becomes dominant, producing 
a bounce to a new Kasner phase and the breakdown of the scheme. These results are 
consistent with the BKL picture for the approach to a generic spacelike singularity. On 
the other hand, for an expanding fluid, our results show that anisotropics die away and 
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Table 2. Asymptotic behaviour of the first-order solution in the particular case 
characterized by the relations H52|) . 



Quantity Behaviour as t — > 



-qT 




H_ 




Hi 
W 


n1iyf')T^-P' hiT 


S3. 

02 


n^{yt')T^-P^ Inr 


0? 


n'iiyf')T^-P-' Inr 



the local evolution is more and more Robert son- Walker-like. This confirms previous 
analyses PIU El E| and is also in line with recent investigations on the evolution of 
non-linearities on superhorizon scales [211 ■ 

The LWAS is based on the assumption that spatial gradients can be neglected with 
respect to time derivatives, whose scale is given by the Hubble parameter. Although 
the approximation scheme is physically well-motivated it turns out that there is not a 
well-established mathematical formulation of the principles on which it is based. The 
difficulty in having such a formulation lies in the difficulty to establish a meaningful way 
of smoothing a clumpy cosmological models, which in turn is related to the averaging 
problem in cosmology, a problem not yet clarified. 

Most of the times the Long- Wavelength Approximation has been directly applied to 
the metric components, and then the spatial gradients of the metric has been neglected, 
which of course is a coordinate dependent procedure. In this work we have applied it to 
the quantities of the covariant fluid approach, which are intrinsically defined once the 
fluid velocity is prescribed, but basically the way in which we have implemented it is 
equivalent to the previous one. The reason for this is that we have neglected completely 
the spatial curvature tensor ^Rab, which means that we have not only neglected spatial 
derivatives of the connection components, Q^F^^, but also the product of them, F^^F^^. 
The second approximation basically means that we have neglected the spatial gradients 
of the metric. 

However, in principle we can generalize this approximation scheme by applying the 
idea at the level of the connection instead of applying it at the level of the metric. That 
is, to neglect all the spatial gradients of the connection, which implies that we will have 
a non-vanishing spatial curvature for the flrst iteration of the scheme. As one can see, 
the equations for this case will have essentially the same form as those of the Bianchi 
models, which are homogeneous models with homogeneous spatial curvature. So we will 
be still inside the spirit of the approximation. The advantage that this generalization 
can have is that it would allow us to include the possibility of having recollapsing regions 
(regions in which > 0), and in this way one could be able to explore a more realistic 
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universe model by using the LWAS, for example studying the evolution, within this 
approach, of a model corresponding to a Robertson- Walker universe with initially small 
perturbations. Since the integration functions of the equations of the LWAS would be 
space dependent this means that we could have regions in which the spatial curvature 
would lead to recoUapse, and other regions that would expand forever. This would give a 
picture closer to the idea we have of the universe from observations, where we have void 
regions that are in expansion and others that are accreting matter and collapsing to form 
cosmic structures. From a mathematical point of view we expect that the generalized 
approach outlined above should be related to the improved LWAS introduced in 
and, for the covariant quantities, should produce a system of ODEs able to describe the 
evolution from expansion to recollapse, much like the equations for "silent" models. 
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